#include <stdlib.h>
#include <errno.h>

#include <assert.h>
#include <petscksp.h>
#include <petscmat.h>

int main(int argc,
	 char** argv)
{
    PetscInitialize(&argc, &argv, (char *)0, (char *)0);
    PetscPrintf(PETSC_COMM_WORLD, "%s", "START\n");

    Vec xVec; /* approx solution, RHS, exact solution */
    Vec xtrueVec;
    Vec bVec;
    Mat matA; /* linear system matrix */
    KSP ksp; /* linear solver context */
    PC pc;  /* preconditioner context */
    PetscReal norm;
    PetscInt its;
    PetscErrorCode ierr;

    /**
     * Reads CRS format into PETSC SEQAIJ format
     */
    PetscInt nr;
    PetscInt nv;

    /* read case0.dat to matrix dimensions */
    FILE* finCaseDat0;
    finCaseDat0 = fopen("case.dat0", "r");
    fscanf(finCaseDat0, "%d", &nr);
    fscanf(finCaseDat0, "%d", &nv);
    assert(nr == 18656);
    assert(nv == 5);
    close(finCaseDat0);
    PetscInt n = nv*nr;
    
    CHKERRQ(ierr = MatCreate(PETSC_COMM_WORLD, &matA));
    
    CHKERRQ(ierr = MatSetSizes(matA,
			       PETSC_DECIDE,
			       PETSC_DECIDE,
			       n,
			       n));
    
    CHKERRQ(ierr = MatSetFromOptions(matA));

    /* Create right-hand-side vector b */
    CHKERRQ(ierr = VecCreate(PETSC_COMM_WORLD, 
			     &bVec));

    CHKERRQ(ierr = VecSetSizes(bVec, 
			       PETSC_DECIDE, 
			       n));

    CHKERRQ(ierr = VecSetFromOptions(bVec));

    CHKERRQ(ierr = VecDuplicate(bVec, &xVec));

    // read RHS Vector B
    double* bArr = (double*)calloc(nv*nr, sizeof(double)) ;

    FILE* fpin = fopen("rhs_value.bin","r");
    int i;
    for(i = 0; i < nr; ++i)
    {
	int m;
	for(m = 0; m < nv; ++m) 
	{
	    fread(bArr+i*nv+m, sizeof(double), 1, fpin);	
	}
    }
    close(fpin);
    
    double* xtrueArr = (double*)calloc(nv*nr, sizeof(double));
    fpin = fopen("x_value.bin","r");
    for(i = 0; i < nr; ++i)
    {
	int m;
	for(m = 0; m < nv; ++m) 
	{
	    fread(xtrueArr+i*nv+m, sizeof(double), 1, fpin);	
	}
    }
    close(fpin);

    PetscInt ix[nr*nv];
    for (i = 0; i < nr*nv; ++i)
    {
	ix[i] = i;
    }

    CHKERRQ(ierr = VecSetValues(bVec, nr*nv, ix, bArr, INSERT_VALUES));
    CHKERRQ(ierr = VecAssemblyBegin(bVec));
    CHKERRQ(ierr = VecAssemblyEnd(bVec));

    //CHKERRQ(ierr = VecSetValues(xtrueVec, nr*nv, ix, xtrueArr, INSERT_VALUES));
    //CHKERRQ(ierr = VecAssemblyBegin(xtrueVec));
    //CHKERRQ(ierr = VecAssemblyEnd(xtrueVec));
    
    // Ax = b
    CHKERRQ(ierr = VecDuplicate(bVec, &xVec));
    CHKERRQ(ierr = VecAssemblyBegin(xVec));
    CHKERRQ(ierr = VecAssemblyEnd(xVec));
        
    int* rowArr = (int*)calloc(nr+1, sizeof(int));

    // read row_ptr data
    fpin = fopen("row_ptr.bin","r");
    for(i=0; i < nr+1; ++i)
    {
	fread(rowArr+i, sizeof(int), 1, fpin) ;
    }
    close(fpin);

    printf("row_ptr.bin has been read. \n");

    // read column index data
    int* colArr = (int*)calloc(*(rowArr+nr), sizeof(int)) ;
    fpin=fopen("col_ind.bin","r");
    int icolst= *(rowArr+0) ;
    int icolen= *(rowArr+nr-1+1);
    for(i = icolst; i <= icolen; ++i)
    {
	fread(colArr+i, sizeof(int), 1, fpin);
    }
    close(fpin);

    printf("col_ind.bin file has been read. \n");

    // read matrix A data
    // double* matArr = (double*)calloc(nv*nv*(*(rowArr+nr-1+1)), sizeof(double));
    
    double* aMat = (double*)calloc(nv*nr*nv*nr, sizeof(double));
    printf("%s\n", strerror(errno));
    fpin=fopen("matrix_coef.bin","r");
    for(i = 0; i < nr*nv; i+=nv)
    {
	icolst = *(rowArr+i);
	icolen = *(rowArr+i+1);
	int j;
	for(j = icolst; j < icolen; ++j)
	{
	    int col = *(colArr+j);
	    int m;
	    for(m = 0; m < nv; ++m) 
	    {
		int n;
		for(n = 0; n < nv; ++n) 
		{
		    fread(aMat+(i*5+m)*nr*nv+col*nv+n, sizeof(double), 1, fpin);
		}
	    }
	}
    }
    close(fpin);

    PetscInt l = nr*nv;
    PetscInt lidx[nr*nv];
    for (i = 0; i < nr*nv; ++i)
    {
	lidx[i] = i;
    }
    CHKERRQ(ierr = MatSetValues(matA, l, lidx, l, lidx, aMat, INSERT_VALUES));

    printf("The matrix A has been read. \n") ;

    CHKERRQ(ierr = MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY));
    CHKERRQ(ierr = MatAssemblyEnd(matA,MAT_FINAL_ASSEMBLY));
  
    CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "%s", "DONE\n"));

    return 0;
    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
       Create the linear solver and set various options
       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    /* 
       Create linear solver context
    */
    CHKERRQ(ierr = KSPCreate(PETSC_COMM_WORLD,
			 &ksp));

    /* 
       Set operators. Here the matrix that defines the linear system
       also serves as the preconditioning matrix.
    */
    CHKERRQ(ierr = KSPSetOperators(ksp,
				   matA,
				   matA,
				   DIFFERENT_NONZERO_PATTERN));
    
    /* 
       Set linear solver defaults for this problem (optional).
       - By extracting the KSP and PC contexts from the KSP context,
       we can then directly call any KSP and PC routines to set
       various options.
       - The following four statements are optional; all of these
       parameters could alternatively be specified at runtime via
       KSPSetFromOptions();
    */
    CHKERRQ(ierr = KSPGetPC(ksp, &pc));
   
    CHKERRQ(ierr = PCSetType(pc, PCJACOBI));
    
    //    ierr = PCFactorSetMatSolverPackage(pc,
    //				       "mumps");
    //CHKERRQ(ierr);

    CHKERRQ(ierr = KSPSetTolerances(ksp,
				    1.e-7,
				    PETSC_DEFAULT,
				    PETSC_DEFAULT,
				    PETSC_DEFAULT));
    
    
    /* 
       Set runtime options, e.g.,
       -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
       These options will override those specified above as long as
       KSPSetFromOptions() is called _after_ any other customization
       routines.
    */
    CHKERRQ(ierr = KSPSetFromOptions(ksp));
    
    /* KSP Set Type to KSPGMRES */
    CHKERRQ(ierr = KSPSetType(ksp, KSPGMRES));
    
    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
       Solve the linear system
       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    /* 
       Solve linear system
    */
    CHKERRQ(ierr = KSPSolve(ksp,bVec,xVec));
     

/*
    VecAssemblyBegin(x);
    CHKERRQ(ierr);
    VecAssemblyEnd(x);
    CHKERRQ(ierr);
*/

/*
    CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "%s", "\n"));
    CHKERRQ(ierr = VecView(x, PETSC_VIEWER_STDOUT_SELF));
    
    CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "%s", "\n"));
    CHKERRQ(ierr = VecView(u, PETSC_VIEWER_STDOUT_SELF));

    CHKERRQ(ierr = PetscPrintf(PETSC_COMM_WORLD, "%s", "\n"));
*/
    /* 
       View solver info; we could instead use the option -ksp_view to
       print this info to the screen at the conclusion of KSPSolve().
    */
    // CHKERRQ(ierr = KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));

    /* MatDestroy(matA) */
    /* VecDestroy(b) */
    /* VecDestroy(x) */

    CHKERRQ(KSPGetIterationNumber(ksp,&its));
    CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,
			"Norm of error %A, Iterations %D\n",
			norm,its));

    PetscFinalize();

    return 0;
}
